samples_surface <- samples %>% filter(surface_sample >= 0)
sample_number <- nrow(samples_surface)
glue::glue("Number of samples: {sample_number}")

otu_long <- otu %>% select(ASV_code, genus, species_wang, species_ASV, contains("OSD")) %>% 
    gather(sample, n_seq, contains("OSD")) %>% filter(sample %in% samples_surface$sample)
otu_long$n_seq <- otu_long$n_seq %>% replace_na(0)

otu_seq <- otu_long %>% group_by(ASV_code) %>% summarise(n_seq_otu = sum(n_seq)) %>% 
    filter(n_seq_otu > 0)

glue::glue("Number of OTUs for all samples: {nrow(otu_seq)}")

glue::glue("Number of reads for all samples: {sum(otu_seq$n_seq_otu)}")

# write_tsv(species_seq, 'OSD Mamiello sequences.tsv')

sample_seq <- otu_long %>% group_by(sample) %>% summarise(n_seq_sample = sum(n_seq))

sample_number_above100 <- nrow(filter(sample_seq, n_seq_sample >= 100))
glue::glue("Number of samples with more than 100 Mamiellophyceae reads: {sample_number_above100}")

otu_long <- otu_long %>% left_join(sample_seq)

summarise_species <- function(otu_long) {
    
    sample_number <- nrow(distinct(as.tibble(otu_long$sample)))
    print(glue::glue("Number of samples taken into account : {sample_number}"))
    
    species_seq <- otu_long %>% group_by(genus, species) %>% summarise(n_seq_species = sum(n_seq))
    
    
    species_samples <- otu_long %>% group_by(species, sample, n_seq_sample) %>% 
        summarise(n_seq = sum(n_seq)) %>% mutate(pct_seq = n_seq/n_seq_sample * 
        100)
    
    species_samples_metadata <- species_samples %>% left_join(select(samples, 
        sample, OSD_station, sample_label)) %>% left_join(metadata) %>% mutate(Latitude = as.numeric(Latitude), 
        Longitude = as.numeric(Longitude))
    
    # Summaries at species level for all species
    
    species_pct <- species_samples %>% group_by(species) %>% summarise(mean_pct_seq_species = mean(pct_seq), 
        max_pct_seq_species = max(pct_seq))
    
    species_presence <- species_samples %>% filter(n_seq > 0) %>% group_by(species) %>% 
        summarise(n_samples_present = n(), pct_samples_present = 100 * n()/sample_number)
    
    species_more_1pct <- species_samples %>% filter(pct_seq > 1) %>% group_by(species) %>% 
        summarise(n_samples_more_1pct = n(), pct_samples_more_1pct = 100 * n()/sample_number)
    
    species_summary <- species_seq %>% left_join(species_pct) %>% left_join(species_presence) %>% 
        left_join(species_more_1pct) %>% filter(!is.na(n_samples_present))  # This to remove the species not present among the selected OTUs
    
    species_stats <- list(summary = species_summary, metadata = species_samples_metadata)
}

otu_long1 <- otu_long %>% rename(species = species_wang) %>% select(-species_ASV)
species_stats_1 <- summarise_species(otu_long1)

species_summary_1 <- species_stats_1[["summary"]]

otu_long1 <- otu_long %>% filter(n_seq_sample >= 100 & !is.na(species_ASV)) %>% 
    rename(species = species_ASV) %>% select(-species_wang)
species_stats_2 <- summarise_species(otu_long1)

species_summary_2 <- species_stats_2[["summary"]]

plot_species_summary <- function(species_summary, file_suffix = "1") {
    
    theme_chlorophyta_map <- theme(legend.position = c(0.1, 0.2), legend.background = element_rect(fill = "transparent"), 
        legend.text = element_text(size = 7), legend.title = element_text(size = 9), 
        plot.tag.position = "topright", plot.tag = element_text(size = 24, face = "bold"), 
        plot.title = element_text(size = 14))
    # Number of sequences per species
    g <- ggplot(species_summary, aes(x = reorder(species, n_seq_species), y = n_seq_species)) + 
        geom_col() + coord_flip() + xlab("Species") + ylab("Total of sequences - all samples")
    print(g)
    
    g_treemap <- ggplot(species_summary, aes(area = n_seq_species, fill = genus, 
        label = species, subgroup = genus)) + ggtitle("OSD - number of sequences") + 
        scale_fill_discrete(name = "Genus") + geom_treemap() + geom_treemap_subgroup_border() + 
        geom_treemap_text(colour = "white", place = "topleft", reflow = T, padding.x = grid::unit(3, 
            "mm"), padding.y = grid::unit(3, "mm"))
    print(g_treemap)
    
    ggsave(plot = g_treemap, filename = str_c("pdf/Treemap_sequences_", file_suffix, 
        ".pdf"), width = 10, height = 7, scale = 2, units = "cm", useDingbats = FALSE)
    
    
    treemap_dv(species_summary, c("genus", "species"), "n_seq_species", "OSD - number of sequences")
    
    
    # Mean pct of sequences per species
    
    g <- ggplot(species_summary, aes(area = mean_pct_seq_species, fill = genus, 
        label = species, subgroup = genus)) + ggtitle("OSD - mean of contribution", 
        subtitle = "") + geom_treemap() + geom_treemap_subgroup_border() + geom_treemap_text(colour = "white", 
        place = "topleft", reflow = T)
    print(g)
    
    treemap_dv(species_summary, c("genus", "species"), "mean_pct_seq_species", 
        "OSD - mean of contribution")
    
    # Graph
    
    graph_samples_present <- ggplot(species_summary, aes(x = reorder(species, 
        pct_samples_present), y = pct_samples_present)) + geom_col() + geom_text(aes(label = n_samples_present), 
        hjust = -0.25) + coord_flip() + xlab("Species") + ylab("% of samples where species detected") + 
        ylim(0, 100) + theme_dviz_grid() + theme_chlorophyta_map + labs(title = "", 
        tag = "A")
    
    
    print(graph_samples_present)
    
    graph_samples_more_1pct <- ggplot(species_summary, aes(x = reorder(species, 
        pct_samples_present), y = pct_samples_more_1pct)) + geom_col() + geom_text(aes(label = n_samples_more_1pct), 
        hjust = -0.25) + coord_flip() + xlab("Species") + ylab("% of samples where species represesents \n more than 1% of Mamiellophyceae") + 
        ylim(0, 100) + theme_dviz_grid() + theme_chlorophyta_map + labs(title = "", 
        tag = "B")
    
    print(graph_samples_more_1pct)
    
    grid_graphs <- gridExtra::grid.arrange(grobs = list(graph_samples_present, 
        graph_samples_more_1pct), ncol = 1, nrow = 2, clip = FALSE, padding = unit(0, 
        "line"))
    
    ggsave(plot = grid_graphs, filename = str_c("pdf/Presence_", file_suffix, 
        ".pdf"), width = 15, height = 20, scale = 1.8, units = "cm", useDingbats = FALSE)
    
    # Relation between samples present and samples more than 1 pct
    g <- ggplot(species_summary, aes(x = pct_samples_more_1pct, y = pct_samples_present, 
        label = species)) + geom_point() + geom_text(size = 2, check_overlap = TRUE, 
        vjust = 1) + xlim(0, 100)
    print(g)
    
    # Relation between sequence per species and samples present
    g <- ggplot(species_summary, aes(x = mean_pct_seq_species, y = pct_samples_present, 
        label = species)) + geom_point() + geom_text(size = 2, check_overlap = TRUE, 
        vjust = 1)
    print(g)
}

plot_species_summary(species_stats_1[["summary"]], file_suffix = "all_ASV")

plot_species_summary(species_stats_2[["summary"]], file_suffix = "selected_ASV")